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The  nonlinear  progressive  wave  equation  (NPE)  is  used  to  investigate  propagation  of  acoustic 
pulses  in  a  shallow  ocean  waveguide.  The  nonlinearity  is  shown  to  affect  transmission  and 
reflection  at  a  fluid-fluid  interface.  It  is  shown  that  one  effect  of  the  nonlinearity  is  to  reduce 
the  critical  angle  for  total  internal  reflection  from  that  of  the  linear  case.  When  loss  versus 
range  for  a  simple  isovelocity  waveguide  is  compared  with  linear  normal  mode  and  parabolic 
equation  calculations,  the  nonlinear  result  shows  increased  transmission  loss  as  long  as  the 
wave  amplitude  is  significant  (an  expected  result  of  shock  processes),  and  increased  bottom 
penetration.  Nonlinear  aging  of  the  waveform  alters  the  excitation  of  linear  modes  in  the 
farfield. 

PACS  numbers:  43.25.Cb,  43.25.Jh,  43.30.Qd 


INTRODUCTION 

The  nonlinear  progressive  wave  equation  (NPE)  mod¬ 
el1  has  recently  been  developed  to  investigate  nonlinear 
acoustic  effects  (including  shocks).  The  model  is  cast  in  the 
time  domain  as  the  most  economical  method  for  represent¬ 
ing  local  nonlinear  phenomena.  For  investigating  propaga¬ 
tion  in  an  ocean  waveguide,  this  model  provides  an  alterna¬ 
tive  to  the  parabolic  equation  ( PE )  and  normal  mode  ( NM ) 
approaches. 

The  NPE  is  derived  from  the  fluid  equations  of  momen¬ 
tum  and  mass  continuity  retaining  lowest-order  nonlinearity 
and  assuming  propagation  within  a  narrow  angle.  The  mod¬ 
el  follows  a  pulse  disturbance  in  a  reference  frame  that  moves 
at  a  constant  representative  sound  speed  c0  in  the  range  di¬ 
rection  (jc) .  The  local  linear  sound  speed  in  the  medium  is 
c  =  c0  +  c,  (x,yj)  where  ct  is  a  small  local  environmental 
departure  from  c0.  Consistent  with  the  PE,  the  NPE  can  be 
derived  from  an  ordering  scheme  in  which  smallness  of  var¬ 
ious  items  is  formally  stated  through  a  scaling  variable  e.  In 
addition  to  the  PE’s  list  of  small  quantities,  the  NPE  scales 
the  wave  amplitude  by  e.  The  list  of  scaled  items  precedes 
Eq.  (2).  The  NPE'  is  stated  as 

— =  -  J-Lr+^r2) 

Dt  dx\  2  ) 

-?±(X  V\Rdx  +  0(e'),  (1) 

2  Jxr 

where  R  —  p'/p0  with  p'  the  acoustic  density  fluctuation  and 
p0  is  the  ambient  density  (assumed  constant  here)  in  the 
medium.  The  coefficient  of  nonlinearity  is  /?=  1  + 
d  log  c/d  log  p  3.5  for  the  ocean.  The  lower  limit  xf  of  inte¬ 
gration  in  ( 1 )  is  located  in  the  quiescent  medium  ahead  of 
the  wave.  On  the  right-hand  side  of  Eq.  (1),  the  c,  term 
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represents  refraction,  the  quadratic  term  contains  nonlinear 
steepening,  and  the  integral  term  includes  diffraction  and 
geometric  spreading. 

The  error  term  in  ( 1 )  is  third-order  in  the  smallness 
parameter  e,  which  scales  each  of  {R,c,,Vj  ,D/Dt },  where 
the  last  two  items  of  the  list  are  the  transverse  Laplacian  and 
the  wave’s  time  derivative  in  the  wave-tracking  frame: 


D  _  d  _d_ 
Dt  dt  +C°  dx' 


(2) 


The  scaling  of  D  /Dt  in  this  context  is  not  an  arbitrary 
assignment,  but  a  fundamental  result  from  the  equations  of 
motion.1  This  is  illustrated  by  the  limit  e-»0,  which  brings 
{c„R,7j }  to  zero  on  the  right  side  of  (1 ).  The  appropriate 
progressive  wave  solution  is  a  linear  plane  wave  in  a  homoge¬ 
neous  medium, /(x  —  c0t),  with  /an  arbitrary  function.  For 
any  solution  of  this  form,  D/Dt  vanishes,  consistent  with  it 
being  scaled  by  e. 

The  equivalence  of  the  NPE  and  PE  under  appropriate 
circumstances  has  been  demonstrated  elsewhere. 1,2  As  with 
the  PE,  any  term  in  ( 1 )  contributes  meaningfully  to  wave 
evolution  only  as  long  as  it  exceeds  the  dominant  error  term. 

One  can  reformulate  the  NPE  in  terms  of  a  dimension¬ 
less  pressure  variable  Q=p7p{pl: 


One  recognizes  (3)  as  differing  from  (1)  only  in  the 
substitution  of  Q  for  R.  The  reason  for  this  symmetry  is  that 
Q—  R  +  Oie2).  Substitution  of  this  expression  into  (3) 
leads  to  an  error  of  order  c’  in  each  term,  consistent  with  ( 1 ) . 
The  pressure  formulation  of  the  NPE  ( 3 )  facilitates  the  han¬ 
dling  of  interfacial  boundary  conditions. 
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In  the  following  sections,  we  present  solutions  of  the 
NPE  for  a  pulse  propagating  along  an  acoustic  waveguide 
descriptive  of  a  shallow  ocean.  We  will  show  how  nonlinear¬ 
ity  affects  transmission  and  reflection  at  the  bottom  (repre¬ 
sented  by  a  fluid-fluid  interface)  and  dispersion  in  the  wave¬ 
guide.  Viscosity  and  sedimentary  attenuation  in  the  ocean 
and  bottom  are  omitted  for  simplicity  (there  is,  however,  a 
wave-absorbing  buffer  zone  below  to  prevent  artificial  grid 
reflections).  The  physical  effects  to  be  investigated  here  are 
not  substantially  altered  by  either  of  these  dissipative  mecha¬ 
nisms.  Linear  results  from  our  model  (0  =  0)  compare  well 
with  existing  linear  results.  Effects  evident  when  nonlinear¬ 
ity  is  included  are  physically  understandable.  Thus  two  nec¬ 
essary  conditions  aie  fulfilled  to  argue  that  the  NPE  is  suit¬ 
able  for  more  complicated  investigations  as  well. 

I.  SIMULATIONS  OF  REFLECTION  AND  TRANSMISSION 
AT  A  FLUID-FLUID  INTERFACE 

In  this  section,  we  present  results  of  simulations  using 
the  pressure  interpretation  of  the  NPE  (3).  We  consider  a 
fluid-fluid  description’  of  the  ocean-bottom  interface  as  in 
Fig.  1.  The  linear  sound  speed  of  the  water  column  is  c0  and 
that  of  the  flat  bottom  is  cb  =  c0  +  c{,  with  c,  >0.  Results 
are  calculated  in  a  reference  frame  moving  at  the  water 
sound  speed  c(),  beginning  with  conditions  illustrated  in  Fig. 
2  ( more  details  below ) .  The  subbottom  has  the  same  proper¬ 
ties  as  the  bottom,  plus  wave-absorbing  numerical  attenu¬ 
ation.  Azimuthal  symmetry  is  assumed. 

Before  giving  the  numerical  results,  we  will  present  a 
simple  theoretical  argument  indicating  how  nonlinearity 
should  affect  total  internal  reflection  at  the  ocean-bottom 
interface.  From  simple  kinematics,  we  can  derive  an  expres¬ 
sion  for  the  critical  grazing  angle  ac  below  which  an  incident 
nonlinear  plane  wave  in  the  ocean  fails  to  excite  a  nonevanes- 
cent  wave  in  the  bottom.  At  the  critical  grazing  angle,  the 
incoming  wave  barely  excites  a  horizontally  propagating 
wave  deep  in  the  bottom.  Without  having  to  solve  for  nonlin¬ 
ear  details  at  the  interface,  a  necessary  condition  for  excita¬ 
tion  of  the  wave  deep  in  the  bottom  is  matching  of  horizontal 
phase  speeds  (i.e.,  a  statement  of  Snell’s  law).  For  this 
matching,  we  need  a  steady-state  incoming  plane  wave.  The 
only  such  wave  with  small  but  finite  amplitude  is  a  weak 
shock  discontinuity.  Mathematically,  this  wave  consists  of 
two  constant  states  Q  =  0  ahead  of  the  shock,  and 
Q  =  Q,  >  0  behind  it.  Within  the  compressed  fluid  behind 
the  shock,  sound  waves  propagate  at  an  increased  speed 

c„,  =c„(  l+/?<2,  )+0(e:),  (4) 


FIG.  I .  Plane-wave  incident  from  ocean  to  bottom,  with  c„  <  ch .  Here,  c„  is 
the  sound  speed  in  the  water,  ch  is  the  sound  speed  in  the  bottom,  or,,,,  is  the 
grazing  angle  of  incidence  (reflection),  and  ar,.r  is  the  angle  of  refraction. 


400  m 
480  m 

600  m 

2  m  202  m 

Range  2  m 

FIG.  2.  Simulation  window  with  initial  waveform.  The  ocean  and  bottom 
are  fluids  of  equal  density  with  sound  speeds  of  1 500  and  1 600  m/s,  respec¬ 
tively.  The  subbottom  has  the  same  properties  as  the  bottom,  plus  numerical 
damping  to  absorb  waves  reflected  from  the  bottom  of  the  numerical  grid. 
The  width  of  the  simulation  window  is  200  m. 


where  subscript  nl  indicates  nonlinear  modification.  The 
shock  front  propagation  speed  (normal  to  the  front )  is  given 
by  the  Rankine-Hugoniot  relation  as  the  average  of  small 
signal  speeds  on  both  sides  of  the  shock,  i.e., 

»,=ebd+l«2.)+0(«2)-  (5) 

The  horizontal  phase  speed  of  this  wave  is  vs  sec  a, .  For 
this  nonlinear  problem,  we  define  a  critical  angle  ac  such 
that  the  incident  wave  excites  a  disturbance  in  the  bottom, 
which,  with  increasing  depth,  approaches  a  horizontally 
propagating  linear  wave  with  phase  speed  c„  +  c,.  Waves 
with  smaller  phase  speeds  are  evanescent  in  the  bottom.  The 
critical  grazing  angle  is  thus  given  by 

c„  +  c,  =  c„  sec  a,.  ( 1  +  \0QX ).  (6) 

For  small  Qs  and  c,/c,„  (6)  gives 

a?s*2(c,/c0) -0QX.  (7) 

This  indicates  a  decrease  of  critical  angle  with  the  ampli¬ 
tude  of  the  incident  wave.  Thus,  a  nonlinear  wave  might  be 
expected  to  disturb  the  bottom  more  than  linear  theory 
would  predict.  In  agreement  with  this,  note  that  for  Qx  >  0 
and  c,  <0coQx/ 2,  there  is  no  real  ac  (no  cutoff  angle  for  the 
nonlinear  wave).  The  sound-speed  discontinuity  is  weak 
enough  that  the  nonlinear  wave  passes  through  to  great 
depth  for  all  grazing  angles,  while  a  linear  wave  would  be  cut 
off  at  a^2(c,/c0). 

Boundary  conditions  for  the  numerical  simulations  are 
as  follows.  The  ocean  surface  is  taken  to  be  a  flat  pressure 
release  surface.  The  lower  boundary  of  the  numerical  simu¬ 
lation  grid  (600-m  depth  in  Fig.  2)  is  a  hard  wall  dt  —  0, 
with  reflected  waves  suppressed  by  a  thick  subbottom  con¬ 
taining  numerical  attenuation  to  be  described  below.  Condi¬ 
tions  that  must  be  satisfied  at  the  ocean-bottom  interface 
(400-m  depth  in  Fig.  2)  are  continuity  of  fluid  pressure  and 
normal  particle  speed.  In  terms  of  the  pressure  variable  Q 
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<?„(*) 

• 

FIG.  3.  The  fluid-fluid  ocean-bottom  inter- 
•  o  face  is  taken  to  be  midway  between  two  verti- 

-  cal  grid  points.  Auxiliary  values  for  finite  dif- 

o  •  ferences  at  the  interface  (open  circles)  are 

generated  from  the  interface  continuity  condi- 
#  tions  (9). 


<?»(*) 


( normalized  on  both  sides  of  the  interface  by  the  bulk  modu¬ 
lus  of  water,  ptpl),  the  conditions  are  continuity  of  Q  and 
p~'d„Q,  where  dn  is  the  derivative  along  the  upward  normal 
to  the  surface. 

One  finite  amplitude  effect  to  be  included  is  deformation 
of  the  interface  by  the  wave.  To  lowest  order  the  slope  of  the 
interface  deformation  £  is 

=  —  f  dzQdx  +  (He2). 

Jx, 

The  deformation  of  the  surface  normal  contributes  a  nonlin¬ 
ear  term  to  the  pressure  continuity  condition  for  the  vertical 
direction: 

PoJdzQu,  =Pol  '9zQb  ~  (Pol'  ~P0h  ')dxQb 

X  ( dzQb  dx  +  0(el),  (8) 

Jx, 

where  subscripts  w  and  b  refer  to  evaluation  on  the  water  and 
bottom  side  of  the  interface,  respectively.  The  fractional  or¬ 
der  in  the  error  term  reflects  that  the  NPE  scales  the  trans¬ 
verse  Laplacian  as  e,  so  that  transverse  derivatives,  e.g.,  dz, 
scale  as  e'n. 

The  interfacial  boundary  conditions  are  represented  nu¬ 
merically  as  follows.  The  undisturbed  interface  is  taken  half¬ 
way  between  two  vertical  grid  points  zs  and  zj+  ,  (Fig.  3). 
The  vertical  grid  may  be  thought  of  as  partitioned  between 
water  and  bottom  variables  Qw  and  Qb .  Artificial  values  Qw 
and  Qh  t  are  generated  from  the  interfacial  boundary  con¬ 
ditions.  In  this  work,  we  examine  effects  of  a  sound-speed 
discontinuity  across  the  interface  without  the  nonlinear 
boundary  term  in  (8).  In  the  results  presented  here,  we  take 
uniform  density  across  the  interface,  p0h  =  p0  ,  so  the  non¬ 
linear  term  in  (8 )  is  automatically  absent.  But  even  allowing 
a  reasonable  density  discontinuity,  neglect  of  the  nonlinear 
term  would  be  justified  by  its  smallness,  estimated  as  follows. 
If p0h  =  l.5p0  ,&ndQ,  is  always  less  than  0.04  (as  will  be  the 
case  considered  below  when  the  wave  first  hits  the  inter¬ 
face),  the  nonlinear  term  on  the  right  side  of  (8)  is  initially 
less  than  1  %  of  the  linear  term,  and  this  ratio  decreases  with 
range.  Unlike  the  nonlinearity  in  (3),  its  effect  does  not  ac¬ 
cumulate  during  propagation.  Retention  of  the  term  would 


affect  results  considerably  less  than  the  uncertainty  in  real 
ocean  bottom  density  values.  Neglecting  the  nonlinear  term, 
the  second-order  finite  difference  statement  of  the  boundary 
conditions  is 

Qwj  +  Qw,,  ,  =  Qb,  +  Qbj  ,  ,  , 

Po~  '(Qbj+ ,  -  Qbj)  =poJ(QwJt ,  -  QWj).  (9) 

Equation  (9)  allows  the  artificial  values  Qw  and  Qb  1  to  be 
eliminated  in  terms  of  the  active  gridpoint  values.  This  rep¬ 
resentation  of  the  interface  conditions  is  physically  consis¬ 
tent  with  the  property  that  when  there  is  no  density  discon¬ 
tinuity  across  the  interface,  the  interface  is  completely 
transparent  (i.e.,  Qw  =  Qb.  and  Qbj+l=  QWj+ ,  )• 

Equation  (3)  is  integrated  in  time  numerically1  using 
Crank-Nicholson  for  the  diffraction  term  and  flux  corrected 
transport  (FCT)  methods4-7  for  the  x  derivative.  The  inter- 
facial  boundary  conditions  are  included  in  the  tridiagonal 
matrix  solution  for  the  Crank-Nicholson  scheme. 

Beneath  the  bottom,  a  damping  zone  is  used  to  absorb 
downward  propagating  waves  that  would  otherwise  reflect 
from  the  bottom  of  the  grid  and  re-enter  the  simulation. 
Damping  of  unwanted  reflections  is  achieved  by  the  addition 
ofaterm  —  Q/rtothe  right-hand  side  of  Eq.  (3),wherethe 
purely  numerical  damping  coefficient  r  varies  with  depth. 

To  study  internal  reflection  at  the  interface,  we  initia¬ 
lized  a  pressure  perturbation  in  the  water  column  and  al¬ 
lowed  it  to  propagate  to  the  bottom.  We  considered  two 
cases:  P  —  0  (linear)  and  p  —  3.5,  where  the  latter  is  an  ap¬ 
proximate  value  for  the  coefficient  of  nonlinearity  in  water. 
The  nonlinear  solutions  we  present  do  not  include  the  effects 
of  cavitation. 

The  initial  configuration  and  parameters  of  the  simula¬ 
tions  are  detailed  as  follows.  The  ambient  fluid  density  is 
taken  to  be  uniform  across  the  interface.  The  initial  pressure 
perturbation  used  in  both  linear  and  nonlinear  simulations  is 
shown  in  Fig.  2.  The  perturbation  0is  a  spherical  wave  with 
maximum  value  Q  =  0. 1  at  a  radius  of  135  m  centered  about 
a  source  point  at  a  depth  of 200  m.  The  wave’s  radial  profile 
is  a  half  sine  wave  of  width  16  m.  A  pressure  perturbation  of 
this  magnitude  ( but  not  of  this  shape )  may  be  expected  from 
an  underwater  explosion.  The  shape  of  the  wave  automati¬ 
cally  relaxes  towards  an  appropriate  self-similar  profile  dur¬ 
ing  propagation.8  The  range  variable  is  x,  and  depth  is  in  the 
negative  z  direction.  The  grid  consists  of  100  by  150  points, 
with  spacing  Sx  —  2  m  and  Sz  =  4  m.  We  assume  a  sound 
speed  c„  =  1 500  m/s  in  the  ocean  and  cb  =  1 600  m/s  in  both 
the  bottom  and  damping  regions.  The  linear  critical  grazing 
angle  (Q,  =  0)  from  (6)  is  20  deg.  The  timestep  is  initially 
0.01  s  and  is  thereafter  adjusted  to  be  no  greater  than  half 
that  allowed  by  the  Courant-Friedrichs-Lewy  (CFL)  sta¬ 
bility  condition.9 

Some  discussion  of  the  starting  field  of  Fig.  2  is  in  order. 
The  NPE,  like  the  PE,  is  a  small-angle  approximation;  yet 
the  starting  field  contains  large-angle  components  with  re¬ 
spect  to  the  range  direction.  The  feature  of  the  NPE  that 
allows  accurate  recovery  of  farfield  results  is  that  the  large- 
angle  wave  components  have  group  velocities  less  than  that 
of  the  advancing  simulation  window  (less  by  a  factor  of  the 
cosine  of  the  propagation  angle).  These  components  pro- 
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gressively  fall  behind  the  field  of  view  and  do  not  affect  the 
farfield  results. 

Simulations  results  are  shown  in  color  in  Fig.  4,  begin¬ 
ning  when  the  maximum  propagation  angle  has  dropped  to 
25  deg  (vertical  scale  compression  by  a  factor  of  3  appears  to 
increase  this  angle).  Pressure  perturbations  in  linear  [4(a)- 


(c)]  and  nonlinear  [4(d)— (f)  ]  cases  are  pictured  at  fixed 
times  for  comparable  values  of  range  ( defined  as  the  distance 
from  the  source  point  to  the  left  of  the  simulation  grid).  The 
range  values  do  not  match  exactly  because  the  timesteps  are 
automatically  adjusted  during  the  runs.  The  color  scale  from 
blue  to  red  differs  from  plot  to  plot.  Reddish  hues  (yellow 


FIG.  4.  Color  contour  plots  of  propagating  pulse  at  fixed  times  in  pulse  following  window.  Yellow  through  red  are  positive,  blue-green  through  violet  are 
negative,  and  green  is  neutral  (zero).  Linear  solution  is  given  in  (a)-(c);  nonlinear  solution  in  (d)-(f).  Distances  from  (he  source  point  to  the  left  of  each 
plot  are:  (a)  452  m.  (b)  1352  m,  (c)  2252  m,  (d)  430  m.  (e)  1338  m,  and  (f)  2207  m. 
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through  red)  are  always  positive,  bluish  hues  (blue-green 
through  violet)  are  negative,  and  green  is  zero.  The  white 
horizontal  line  through  the  simulation  region  is  the  bottom. 
As  the  wave  propagates  far  from  the  source,  the  effective 
grazing  angle  with  the  bottom  becomes  approximately  in¬ 
versely  proportional  to  range.  In  this  regard,  our  results  at 
increasing  range  show  effects  of  decreasing  grazing  angle. 

Figure  4(a)  and  (d)  is  constructed  from  simulation 
data  at  452  and  430  m,  respectively,  with  (a)  being  linear 
and  (d )  nonlinear.  The  linear  wave  intersects  the  bottom  at  a 
grazing  angle  of  approximately  25  deg.  Since  the  linear  criti¬ 
cal  grazing  angle  is  20  deg,  energy  can  still  propagate  into  the 
bottom  for  a  short  while.  An  obvious  difference  between  the 
two  cases  is  the  shape  of  the  waveforms.  The  propagating 
wave  in  the  linear  case  maintains  its  essential  half-sine  wave 
profile,  while  the  nonlinear  wave  has  steepened  into  a  shock 
ahead  of  what  appears  to  be  an  exponential  tail. 

Figure  4(b)  and  (e)  shows  the  pressure  perturbations  at 
1352  and  1338  m,  respectively.  The  linear  wave’s  grazing 
angle  with  the  bottom  (approximately  8  deg)  is  now  well 
below  critical,  and  propagation  into  the  bottom  has  been  cut 
off.  The  nonlinear  wave  continues  to  propagate  as  a  sharp 
front,  while  the  linear  wave  maintains  a  rounded  profile.  At 
this  range,  the  bottom-reflected  linear  wave  is  a  positive 
pulse  followed  by  a  negative  pulse  of  approximately  equal 
amplitude.  The  nonlinear  wave  has  a  predominantly  positive 
reflection.  There  is  a  local  maximum  in  linear  and  nonlinear 
waves  near  the  bottom  evident  as  a  red  area.  Comparing  the 
shape  and  size  of  this  region  reveals  that  much  of  the  nonlin¬ 
ear  wave’s  energy  resides  in  the  bottom-reflected  wave,  while 
the  linear  wave’s  energy  is  primarily  in  the  direct  arrival. 
This  is  a  result  of  the  nonlinear  direct  arrival  being  weakened 
by  shock  processes,  rather  than  a  strengthening  of  the  re¬ 
flected  wave.  (The  color  plot  automatically  rescales  each 
frame  to  minimum  and  maximum  values  for  clarity  of  pre¬ 
sentation.  ) 

Figure  4(b)  and  (e)  shows  the  transient  behavior  of  the 
energy  that  penetrated  the  bottom  earlier  in  the  wave’s  evo¬ 
lution.  This  energy  appears  as  a  finger  somewhat  ahead  of 
the  main  wave.  Separation  from  the  wave  occurs  when  the 
grazing  angle  goes  subcritical.  From  this  point  on,  addi¬ 
tional  energy  does  not  radiate  into  the  bottom,  but  runs 
ahead  and  radiates  back  into  the  water  column  ( this  is  evi¬ 
dent  from  movies  made  in  a  similar  calculation).  Energy 
diagnostics  (not  shown  here)  confirm  that  despite  the  non¬ 
linear  wave’s  weakening  by  shock  processes,  it  has  deposited 
more  energy  into  the  bottom  by  this  time  than  the  linear 
wave.  We  interpret  this  as  being  due  to  a  smaller  critical 
angle  for  the  nonlinear  wave,  and  discussed  above. 

Finally,  we  compare  the  two  cases  at  2252  and  2207  m 
downrange,  as  seen  in  Fig.  4(c)  and  (f),  respectively.  Dis¬ 
turbance  of  the  bottom  by  the  linear  wave  (c)  is  negligible, 
while  the  nonlinear  wave  (f)  still  shows  some  disturbance. 
The  bottom  reflection  of  the  linear  wave  has  shifted  in  phase 
and  become  mostly  negative,  as  predicted  by  linear  theory. 
This  leads  to  destructive  interference  at  and  below  the  inter¬ 
face.  The  nonlinear  reflected  wave,  however,  is  still  substan¬ 
tially  positive,  resulting  in  constructive  interference  with  the 
direct  arrival  (this  mechanism  for  increased  bottom  distur¬ 


bance  is  separate  from  the  nonlinear  decrease  in  the  critical 
angle).  Another  physical  effect  visible  in  Fig.  4  is  that  the 
nonlinear  wave  propagates  slightly  faster  than  the  linear 
wave.  By  the  time  of  Fig.  4(c)  and  (f),  the  forward  edge  of 
the  nonlinear  wave  has  crept  forward  relative  to  the  coordi¬ 
nate  system  moving  at  speed  c,„  while  that  of  the  linear  wave 
has  remained  fixed. 

The  above  comparison  can  be  summarized  as  follows.  A 
linear  compressional  wave  incident  on  a  fluid-fluid  interface 
undergoes  a  transition  from  positive  to  negative  reflection  as 
the  effective  grazing  angle  drops  below  the  critical  value  and 
tends  to  zero.  When  nonlinear  effects  are  included,  the  criti¬ 
cal  angle  is  smaller,  the  reflected  wave  is  primarily  positive, 
the  disturbance  of  the  bottom  is  greater,  and  pulse  propaga¬ 
tion  is  slightly  faster. 

II.  DISPERSION  OF  A  WAVE  PACKET 
A.  Normal  modes 

The  results  of  the  previous  section  suggest  that  large- 
amplitude  waves,  such  as  those  resulting  from  an  under¬ 
water  explosion,  will  interact  differently  with  the  ocean-bot¬ 
tom  interface  than  small-amplitude  waves,  because  of  the 
nonlinear  contribution  to  the  local  sound  speed.  Because  of 
this  effect,  we  expect  the  dispersion  of  a  nonlinear  pulse  to  be 
different  than  in  the  linear  case.  In  this  section,  we  present 
numerical  results  using  the  NPE  model  designed  to  high¬ 
light  some  of  those  differences. 

We  consider  an  idealized  shallow  ocean  like  that  pic¬ 
tured  in  Fig.  2.  In  this  case,  we  assume  an  ocean  depth  of  200 
m  together  with  200  m  of  bottom,  the  lower  four-fifths  of 
which  is  progressively  damped  to  avoid  numerical  reflection 
from  the  lower  simulation  boundary.  At  a  source  depth  of  50 
m,  we  initialize  a  five-cycle  sine-wave  packet  with  a  wave¬ 
length  consistent  with  a  frequency  of  50  Hr  and  limit  the 
angular  spread  to  38  deg.  A  Gaussian  envelope  of  1/e  width 
1 .77  cycles  multiplies  the  sine-wave  packet.  The  sound  speed 
in  the  water  column  is  taken  to  be  1500  m/s  and,  in  the 
bottom,  is  1 550  m/s.  The  initial  perturbation  is  then  allowed 
to  propagate  downrange. 

Data  from  a  simulated  array  of  15  hydrophones  placed 
at  various  downrange  locations  are  stored  during  the  run  as  a 
collection  of  time  series.  We  ran  a  linear  case  (P  =  0)  out  to 
20  km  and  repeated  the  simulation  for  a  nonlinear  case 
(P  =3.5).  An  initial  amplitude  £)max  =  0.2  was  chosen  to 
illustrate  differences  between  linear  and  nonlinear  propaga¬ 
tion. 

Normal  mode  solutions  <i>  to  the  wave  equation  for  a 
harmonic  point  source  of  a  single  frequency  at  depth  z„  (for 
the  case  of  azimuthal  symmetry)  are 

<f>(r,z)  =  ~p(z(t)^?u„  {z{))un(z)H[')(k„r),  (10) 

II 

where p(z(t)  is  the  density  at  the  source  depth  z„  and  H ’  is 
the  Hankel  function  of  the  first  kind.  The  normal  modes 
u„  ( z )  are  eigenfunctions  of  the  equation 

d2u„(z)/dzr  +  [Ar2(z)  -k2„]u„(z)  =0.  (11) 

In  Eq.  (11),  k(z)  =  o/c{z)  is  the  wavenumber  and  k  \ 
is  the  eigen  value  for  mode  n.  The  normal  modes  propagate  at 
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FIG.  5.  Normal  mode  solution  from  SNAP.  The  three  modes  are  shown. 

different  speeds  and  will  separate  as  the  pulse  propagates 
downrange.  The  normal  modes  of  the  waveguide  have  been 
calculated  numerically10  and  results  are  depicted  in  Fig.  S. 

Figure  6  shows  simulated  hydrophone  time-series  plot¬ 
ted  as  a  function  of  reduced  time  using  the  NPE  model  for 
the  linear  case.  Reduced  time  is  defined  for  each  vertical 
array  relative  to  the  time  at  which  the  pulse  first  arrives  at 
the  hydrophone  array.  Figure  6(a)  shows  the  readings  at  a 
range  of  600  m,  and  clearly  there  is  no  mode  separation. 


Figure  6(b)  shows  the  readings  for  the  hydrophone  array  at 
S  km  downrange.  At  this  distance,  the  mode  separation  is 
incomplete,  although  there  is  some  indication  of  modes  two 
and  three  splitting  away  from  the  packet.  Further  down- 
range  at  20  km,  the  modes  have  separated,  as  seen  in  Fig. 
6(c).  The  dashed  lines  in  Fig.  6(c)  indicate  approximate 
arrival  times  of  normal  mode  maxima.  We  can  see  that  the 
general  appearance  of  the  modes  is  consistent  with  the  pre¬ 
dictions  of  normal  mode  analysis.  For  example,  in  the  sec¬ 
ond  mode  arrival  [reduced time 0.20 to 0.25 sin  Fig.  6(c)]  a 
phase  reversal  is  evident  along  the  vertical  dashed  lines,  with 
a  null  at  approximately  108-m  depth.  Figure  6(d)  shows  the 
amplitudes  at  each  of  the  dashed-line  arrival  times  as  a  func¬ 
tion  of  hydrophone  depth. 

The  simulated  hydrophone  data  in  the  nonlinear  case 
reveal  differences  from  the  linear  results.  Figure  7(a)  shows 
these  results  for  a  range  of  600  m.  The  most  apparent  differ¬ 
ence  is  the  wave  steepening.  The  individual  waveforms  are 
also  more  complex.  Even  though  the  presence  of  nonlinear¬ 
ity  does  not  allow  the  rigorous  mathematical  extraction  of 
normal  modes,  we  expect  similarity  between  the  cases  since 
the  nonlinearities  considered  here  are  weak.  After  an  initial 
phase  in  which  the  nonlinear  wave  loses  energy  to  shock 
processes  and  increased  bottom  penetration,  its  interaction 
with  the  waveguide  becomes  essentially  linear. 


RANGE  0.6  km  (a) 


RANGE  5  km  (b) 


- aM/V- — 


-'vyyv\A/'— 
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REDUCED  TIME  T-R/1.8  (S) 

RANGE  20  km  (c) 


REDUCED  TIME  T-R/l.S  (S) 

(d) 


,  FIG.  6.  Numerically  simulated  hydro¬ 
phone  readings  at  ranges  of  (a)  600  m,  (b) 
5  km,  and  (c)  20  km  for  the  linear  case.  In 
(c),  the  nine  vertical  lines  are  used  to 
point  out  the  mode  structure  of  the  simu¬ 
lated  hydrophone  data.  At  a  range  of  20 
km,  there  is  a  clear  separation  of  the 
modes.  In  (d),  the  modes  marked  in  (c) 
are  displayed. 


REDUCED  TIME  T-R/1.8  (8) 
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RANGE  0.6  km 


(a) 


RANGE  5  km  (b) 


REDUCED  TIME  T-R/l.5  (S)  REDUCED  TIME  T-R/1.5  (S) 

FIG.  7.  Nonlinear  case  for  the  same  condi¬ 
tions  as  Fig.  6. 


RANGE  20  km  (c)  (d) 


REDUCED  TIME  T-R/1.5  (S) 


With  these  points  in  mind,  we  examine  the  hydrophone 
data  from  the  simulation  of  the  nonlinear  wave  shown  in  Fig. 
7(c)  at  20  km  downrange.  Although  the  results  differ  from 
Fig.  6(c),  we  can  make  out  the  arrival  of  the  three  “modes” 
as  indicated  by  the  dashed  lines.  Plotting  the  hydrophone 
signal  amplitude  along  the  dashed  lines  yields  Fig.  7(d). 
Again  there  is  relatively  good  agreement  between  these 
modes  and  th«  results  in  Fig.  5.  An  exception  is  mode  1 
identified  in  the  first  three  dashed  lines  of  Fig.  7(c)  and 
plotted  in  (d).  While  the  mode  profile  in  the  linear  case 
agrees  favorably  with  that  of  Fig.  5,  the  nonlinear  profiles  are 
different,  exhibiting  an  amplitude  maximum  at  a  depth  of 
about  80  m  rather  than  the  deeper  1 10  m  predicted  by  nor¬ 
mal  mode  analysis.  Since  these  slices  in  depth  are  taken  near 
the  front  of  the  propagating  signal,  this  departure  from  the 
linear  result  may  be  an  indication  of  higher  early  bottom 
losses  in  the  nonlinear  case  as  compared  to  the  linear  case. 

The  results  of  the  normal  mode  comparison  can  be  sum¬ 
marized  as  follows.  Agreement  is  good  between  the  eigen- 
modes  predicted  by  linear  theory  and  the  separate  arrivals 
that  emerge  during  propagation  of  a  wave  packet  in  the  lin¬ 
ear  NPE  calculation.  In  the  nonlinear  case,  time  series  of  the 
amplitude  sampled  near  the  source  depth  at  a  relatively 
short  distance  downrange  exhibit  the  steepening  expected  of 
a  weak  shock  wave.  As  the  nonlinear  pulse  propagates,  the 
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combined  effects  of  shock  dissipation  and  more  efficient  cou¬ 
pling  with  the  bottom  reduce  the  amplitude  relative  to  the 
linear  wave.  Far  downrange,  dispersion  begins  to  separate 
the  wave  into  normal  modes  recognizable  from  linear  theo¬ 
ry.  The  differences  that  persist  between  linear  and  nonlinear 
cases  are  evidence  of  nonlinear  aging."  (A  nonlinear  wave 
propagating  in  a  nondispersive  medium  undergoes  changes 
in  its  profile  and  in  the  shape  of  its  frequency  spectrum. 
These  changes  are  “remembered”  when  the  wave  weakens 
enough  to  be  considered  linear. ) 

B.  Frequency  spectra 

In  the  previous  section,  we  examined  normal  modes  of 
the  waveguide  by  sampling  time  series  from  vertical  arrays 
of  points  at  several  well-separated  locations  downrange. 
Fourier  analysis  of  time  series  at  the  source  depth  yields  in¬ 
formation  about  changes  in  the  frequency  spectrum  with 
range. 

Figure  8  shows  frequency  spectra  for  the  simulations  of 
Sec.  II A  at  specific  locations  downrange  for  the  linear  case. 
Throughout  the  range  from  600  m  to  20  km,  the  spectrum 
remains  strongly  peaked  near  the  initial  central  frequency  of 
50  Hz.  The  Gaussian  envelope  has  its  peak  at  0  Hz  and  along 
with  the  50-Hz,  five-cycle  sine  wave  forms  harmonics  at  50- 
Hz  intervals  (i.e.,  100, 150,  200  Hz,  etc.). 
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FIG.  8.  Spectrum  of  propagating  pulse  at  ranges  (a)  600  m,  (b)  5  km,  (c) 
10  km,  and  ( d )  20  km.  Ti.e  pulse  used  is  a  five-cycle  sine  wave  multiplied  by 
a  Gaussian  envelope. 


Figure  9,  depicting  the  nonlinear  results,  presents  a 
striking  contrast.  While  the  predominant  frequency  remains 
near  50  Hz,  the  sidelobes  containing  other  frequencies  are  of 
considerably  greater  amplitude  than  in  the  linear  case,  par¬ 
ticularly  at  the  short  range  of  600  m.  This  is  the  stage  at 
which  the  wave  steepens  into  a  shock,  giving  rise  to  high- 
frequency  components.  As  the  wave  continues  to  propagate, 
the  composition  of  frequencies  changes  and  begins  more 
closely  to  resemble  the  linear  case. 

C.  Loss  versus  range 

In  this  section,  we  examine  loss  versus  range  for  a  single 
frequency,  benchmarking  the  NPE  against  standard  linear 


Range  (km) 

FIG.  10.  Loss  versus  range  for  the  linear  simulations,  (a)  The  NPE 
(0  =  0)  (solid)  versus  the  normal-mode  (dashed)  calculations,  while  (b) 
shows  the  PE  (solid)  versus  the  normal-mode  (dashed)  calculations.  Both 
agree  favorably  with  the  normal-mode  solution. 


techniques  in  the  linear  regime.  We  also  compare  the  nonlin¬ 
ear  case  with  linear  predictions. 

The  50-Hz  spectral  component  from  Sec.  II  B  is  plotted 
against  range  in  Fig.  10.  Shown  here  are  the  NPE  versus  the 
normal  mode  solutions  [in  Fig.  10(a)]  and  the  PE  results 
from  the  PAREQ  model 12  versus  the  normal  mode  solutions 
[in  Fig.  10(b)  ].  We  see  that  both  the  NPE  and  PE  models 
are  in  good  agreement  with  the  normal  mode  predictions. 

Figure  1 1  is  a  comparison  of  losses  between  the  linear 
normal  mode  solution  and  the  NPE  results,  including  non¬ 
linearity.  There  is  a  substantial  initial  loss  of  energy  that  we 
attribute  to  shock  processes  and  to  increased  bottom  pene¬ 
tration.  Later,  the  losses  parallel  those  in  the  linear  case. 


Frequency  (Hz) 


Frequency  (Hz) 


FIG.  9.  Nonlinear  case  for  the  same  conditions  as  Fig.  8. 


FIG.  1 1 .  Loss  versus  range  comparing  the  nonlinear  (0  =  3.5,  Qmm%  =  0.2 ) 
NPE  calculation  (solid  line)  with  the  linear  normal-mode  result  (dashed 
line). 
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III.  SUMMARY  AND  CONCLUSION 

We  have  used  the  NPE  model  to  compare  linear  and 
nonlinear  acoustic  propagation  in  a  shallow-water  wave¬ 
guide.  Nonlinear  effects  evident  in  the  comparison  are:  (1 )  a 
smaller  critical  grazing  angle  at  the  bottom,  resulting  in  en¬ 
hanced  transmission  of  wave  energy;  (2)  energy  loss  near 
the  source  attributed  to  shock  formation;  and  (3)  nonlinear 
aging  evident  in  differences  in  normal-mode  excitation  as  the 
wave  weakens.  In  its  linear  mode,  the  NPE  was  compared 
with  parabolic  equation  and  normal-mode  methods.  Trans¬ 
mission  loss  versus  range  from  all  three  are  in  good  agree¬ 
ment. 

Using  simulated  data  from  sample  points  arranged  as 
hydrophone  arrays,  we  have  been  able  to  identify  eigen- 
modes  of  the  waveguide  and  compare  them  with  normal¬ 
mode  results.  In  the  linear  case,  the  simulated  profiles  agree 
with  mode  predictions.  In  the  nonlinear  case,  differences  are 
evident,  although  similarities  with  the  linear  eigenmodes  are 
identifiable.  Of  the  three  nonlinearly  distorted  modes,  the 
first  two  exhibit  some  departure  from  the  calculated  linear 
modes  and  appear  to  peak  nearer  the  surface.  An  example  of 
a  nonlinear  mode  that  agrees  with  the  linear  result  is  also 
evident.  This  is  consistent  with  the  idea  that  the  nonlinear 
pulse  loses  a  large  amount  of  energy  to  shock  processes  and 
to  increased  bottom  penetration  (lower  critical  angle).  After 
the  initial  nonlinear  reduction  in  wave  amplitude,  the  signal 
propagates  and  disperses  in  a  manner  similar  to  the  linear 
case.  Differences  that  persist  after  the  wave  weakens  [  Figs. 
6(c)  and  (d)  vs  7(c)  and  (d)]  are  attributed  to  nonlinear 
aging  of  the  wave. 

Data  from  sample  points  were  also  Fourier  analyzed  to 
obtain  frequency  spectra  during  propagation.  There  are 
strong  differences  shown  in  this  analysis  between  linear  and 
nonlinear  waves.  The  steepening  of  the  shock  front  is  evident 
in  the  generation  of  high-frequency  components.  The  result¬ 
ing  spectrum  downrange  is  much  broader  than  for  the  corre¬ 
sponding  linear  case,  and  the  central  frequency  is  also  lower. 


Finally,  the  consistency  of  our  linear  results  with  the 
existing  literature,  taken  together  with  the  physical  plausi¬ 
bility  of  our  nonlinear  results,  lends  credibility  to  the  NPE  as 
a  suitable  numerical/theoretical  tool  for  studying  nonlinear 
acoustic  phenomena. 
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